#!/usr/bin/env python
# -*- coding: utf-8 -*-
# @Time    : 2021/3/28 14:28
# @Author  : Rem~
# @File    : 004SineSignal_WithWindow.py
# @function: 在003的基础上加窗
# 模拟真实的测量信号

import numpy as np
import matplotlib.pyplot as plt

# sample frequency
fs = 100
# signal frequency
f = 5
# signal phase
fai = 55*np.pi/180
# sample point
N = 1000
# x axis
x = np.linspace(0, N-1, N)
# amplitude
amp = 2**12
bias = amp
# 探测器接收的光信号
sig = amp*np.cos(2*np.pi*f*(x/fs)+fai)+bias
# 高斯窗函数
wn = np.exp(-((x-N/2)/N*4)**2)
# 高斯窗下的光信号
sig_wn = sig*wn

fig = plt.figure("Time&Frequency Domain Analyse")
axes1 = fig.add_subplot(3, 2, 1)
axes2 = fig.add_subplot(3, 2, 2)
axes3 = fig.add_subplot(3, 1, 2)
axes4 = fig.add_subplot(3, 3, 7)
axes5 = fig.add_subplot(3, 3, 8)
axes6 = fig.add_subplot(3, 3, 9)

axes1.plot(x, sig, "-b")
axes1.set(xlabel="sample point", ylabel="Time Domainsignal amplitude(LSB)")
axes1.grid(True)

axes2.plot(wn)
axes2.set(xlabel="sample point", ylabel="Gauss Window amplitude(LSB)")
axes2.grid(True)

axes3.plot(sig_wn)
axes3.set(xlabel="sample point", ylabel="Signal in Gauss Window amplitude(LSB)")
axes3.grid(True)

sig_wn_fft = np.fft.fft(sig_wn)
sig_wn_fft_magnitude = np.abs(sig_wn_fft)/N*2
sig_wn_fft_phase = np.angle(sig_wn_fft)*(180/np.pi)

axes4.stem(x[0:int(N/2)-1]*(fs/N), sig_wn_fft_magnitude[0:int(N/2)-1])
axes4.set(xlabel="frequency axis(MHz)", ylabel="signal amplitude(LSB)")
axes4.grid(True)

axes6.plot(x[0:int(N/2)-1]*(fs/N), sig_wn_fft_phase[0:int(N/2)-1])
axes6.set(xlabel="frequency axis(MHz)", ylabel="signal phase(°)")
axes6.grid(True)

sig_fft = np.fft.fft(sig)
sig_fft_magnitude = np.abs(sig_fft)/N*2
sig_fft_phase = np.angle(sig_fft)*(180/np.pi)

wn_fft = np.fft.fft(wn)
wn_fft_magnitude = np.abs(wn_fft)/N*2
wn_fft_phase = np.angle(wn_fft)*(180/np.pi)
axes5.plot(x[0:int(N/2)-1]*(fs/N), sig_fft_phase[0:int(N/2)-1], 'b')
axes5.plot(x[0:int(N/2)-1]*(fs/N), wn_fft_phase[0:int(N/2)-1], 'r')
axes5.set(xlabel="frequency axis(MHz)", ylabel="signal phase(°)")
axes5.grid(True)

sig_wn_fft_abs_max = np.max(sig_wn_fft_magnitude[5:int(N/2)-1])
sig_wn_fft_abs_max_f = np.argmax(sig_wn_fft_magnitude[5:int(N/2)-1])+5
sig_wn_fft_max_f_phase = sig_wn_fft_phase[sig_wn_fft_abs_max_f]
wn_fft_max_f_phase = wn_fft_phase[sig_wn_fft_abs_max_f]
sig_fft_max_f_phase = sig_fft_phase[sig_wn_fft_abs_max_f]
print('频谱峰值:', sig_wn_fft_abs_max, '\n最大峰值点:', sig_wn_fft_abs_max_f*(fs/N), 'MHz',
      '\n最大频率点相位', sig_wn_fft_max_f_phase, '°')
print('sig最大频率点相位:', sig_fft_max_f_phase, '°\n',
      'wn最大频率点相位:', wn_fft_max_f_phase, '°')

# fig.tight_layout()
# show放print前面时，console界面不会出现相应的print信息
plt.show()
